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- - - We examine the stability of hadron resonance gas models by extending them to take care of 

undiscovered resonances through the Hagedorn formula. We find that the influence of unknown 
resonances on thermodynamics is large but bounded. Hadron resonance gases may be internally 
consistent up to a temperature higher than the cross over temperature in QCD; but by examining 
quark number susceptibilities we find that their region of applicability seems to end substantially 
O . below the QCD cross over. We model the decays of resonances and investigate the ratios of particle 

D ■ yields in heavy-ion collisions. We find that observables such as hydrodynamics and hadron yield 

' ratios change little upon extending the model. As a result, heavy-ion collisions at RHIC and LHC 

are insensitive to a possible exponential rise in the hadronic density of states, thus increasing the 
stability of the predictions of hadron resonance gas models in this context. 



o 
o 

(N 



(N 
(N 



PACS numbers: 



P-i 

I ' Hadron resonance gas models (HRGM) in their modern form were first explored in [1-4], where the emphasis was 

O^' on exploring the phase diagram. Then, following the CERN SPS experiments with heavy-ions, it was found that 

^ integrated particle yields could be explained in HRGMs [5-7] . Typically, one defines the thermodynamics of HRGMs 

' through the summed free energy — 



log Z{V,T,^l)^ / dm[pM{m)logZt,{m,V,T, ^i) + pB{m)\ogZf{m,V,T, ^,)] , (1) 
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> where the gas is contained in a volume V, has temperature T and chemical potential p., Zt is the partition function 
for an ideal gas of bosons with mass m, Zf of an ideal gas of fermions, puiiTi) is the spectral density of mesons, and 



. psim) of baryons. From this one computes the energy density, E, by taking a derivative with respect to 1/T, and the 
■ pressure, P, by taking a derivative with respect to V. One can also find the conformal symmetry breaking measure 
O ; (E- 3P)/T^, the entropy density, S/T^ = {E + P)/T'^, the speed of sound, = dP/dE, and the specific heat, Cv 
Hadron properties enter these models through pb,m and the treatment of the decay width and size of the hadron. 
We follow [7] in neglecting an excluded volume effect for hadrons, i.e., we treat the hadrons as effectively point-like. 
In constructing the thermodynamics, the decay width is usually neglected; we accept this assumption through the 
part of this work which deals with the equation of state, but, as usual, include decays when we explore other aspects 
of HRGMs. One aspect common to models of the type explored in [5-7] is to take the observed spectrum of hadrons 
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d up to some cutoff A, i.e., write 
_ _ I 

mi<A 

PM.,B{ni)= ^ giS{m~mi), (2) 

i 

where are the masses of the known hadrons and gi the degeneracy factor (2 -I- l)(2/i -I- 1) where Ji is the spin and 
li the isospin, and the sum is over meson or baryon states, as appropriate. We call such models HRGl. We include 
in the sum above all the states which are reasonably well established, and labeled as better than 1 star in the particle 
data book. 

Clearly the predictions of thermodynamic quantities in HRGl depend on the mass cutoff A. This is shown in Figure 
1. Even at temperatures as low as 150 MeV, there is a 5% increase in the energy density when A is increased from 
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FIG. 1: Thermodynamic quantities in the hadron resonance gas variant which we caU HRGl. As the mass cutoff, A on the 
states included in the model is changed the predictions change monotonically. These are results at zero chemical potential. 



1600 MeV to 1800 MeV, and a further 2% increase when it is increased to 2000 MeV. Changes at temperatures of 
170 MeV or so are significantly larger. Predictions of quantities such as the speed of sound, c^, and the specific heat, 
Cv, are significantly less stable. Recent quark model [8] and lattice computations [9] lead us to believe that there is 
a much higher density of hadronic states in the mass range 2-3 GeV than below 2 GeV. If so, it is not clear by how 
much the predictions will change as A is increased to 3 GeV or beyond. 

In order to explore the stability of predictions from HRGMs, we develop a variant of these models in which we 
take the observed spectrum of states up to a certain cutoff A, as before, and above this we put in an exponentially 
rising cumulative density of hadron states in [10, 11]. It is well known that such a density of states overcomes the 
exponential suppression of higher mass resonances in eq. (1). As a result, there is a limiting temperature in these 
models, which was, long ago, taken as evidence of a QCD phase transition [12]. So, the extended model, IIRG2, has 

T?li<A 

Ph{m) = V gAm - m,) + ^e'"/^-e(m - A), (3) 

i 

where the model parameters Th and an are fitted to data on the cumulative distribution of different sets of hadrons, 
h. Similar models have been used to study observables as diverse as dilepton production [13] to chemical equilibration 
[14] . A comparison of the equation of state with different density of states for the Hagedorn spectrum is presented 
later. 

We perform fits separately to non-strange mesons, strange mesons, non-strange baryons and strange baryons. All 
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FIG. 2: Fit of the parameters of the Hagedorn spectrum to states from the particle data book for masses up to 2000 MeV. The 
band is wide enough to include all but the two extreme points in the cumulative distribution for each of four classes of hadrons. 



Class (/i) 


mill 


max 


Mesons (non strange) 


0.81 


1.64 


Mesons (strange) 


0.64 


0.94 


Baryons (non strange) 


0.73 


1.47 


Baryons (strange) 


0.50 


1.10 



TABLE I: The extreme values of au such that for Th = 340 MeV the cumulative density of hadron states up to 2000 MeV lies 
within the bands except for the two extreme points on each side. 



states in the particle data book with rating higher than 1 star and mass up to 2000 MeV are used in these fits. The 
quantity Th is obtained by a global fit to all four types of hadrons, whereas we fit upper and lower values for ah for 
each class of hadrons. We have allowed the band of ah/Tn to be wide enough that all but two of the extreme high 
and low points are within the band. The goodness of such fits can be gauged from Figure 2. The fitted parameter 
sets to be used in eq. (3) are given in Table I. Note that the best fit value of Th is significantly higher than Tc, the 
cross over temperature of QCD [15]. This is not unexpected, since there is good evidence that Th and T^ are not 
equal, at least at finite [20]. As a result, the model in eq. (3) is internally consistent well beyond Tc- 

Using this extended model, we get the predictions of thermodynamic quantities shown in Figure 3. Note that for 
each value of A there is an upper and lower limit to the prediction for any thermodynamic quantity, one from the 
maximum and the other from the minimum of the at- As one changes A, the band for higher A lies entirely within that 
for lower A. This pleasant property implies a stability of the predictions of the HRG2 model for bulk thermodynamic 
quantities. Even for and Cv, predictions of HRG2 are stable in this sense up to a temperature of 160 MeV. The 
allowed band at 150 MeV is about 2-7% for A = 2000 MeV. 

A different density of states for the Hagedorn model has also been used in the literature [21, 22]. The density of 
states in the HRG2 with this change would be 

mi<A 

p,(m) = 9^Sim - m,) + ^^^^e"/^«e(m - A). (4) 

Typically this model is used with mp = 500 MeV. We have used this canonical value as well mp = 250 MeV and 1000 
MeV. Note that in the limit as m^O oo, this model reduces to the one in eq. (3). The best fit value of Th changes 
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FIG. 3: Thermodynamic quantities in the hadron resonance gas variant which we call HRG2. At each value of the parameter 
A, there are upper and lower limits on every thermodynamic quantity. The predictions of HRG2 are stable, since the allowed 
band shrinks as A increases, and the band for higher A is entirely contained inside that for a lower A. If QCD has a cross over 
at finite temperature, and the same equation of state is valid on both sides, then falls monotonically across this point, and 
the softest point is not at Tc. 
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FIG. 4: Comparing predictions of thermodynamics using HRG2 with the parametrisation of eq. (3) (full line) and eq. (4) 
(dashed line). The uncertainty bands of all quantities overlap over the range of temperatures of interest. 
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significantly as wc vary mo, being 210 MeV for the central value of too, and sliding up to 250 McV as toq changes 
to 1000 MeV. The quality of fit improves marginally as toq increases, being best for the density of states in eq. (3). 
The uncertainty in the value of the derived hadronic quantity, Th, is closely related to the fact that a string model 
of hadrons is not a unique and self-consistent theory. As a result, the pre-exponential factor can be tweaked at will, 
resulting in large possible changes to the string tension, or, equivalently, to Th- 

Using Too = 500 MeV and Th = 210 MeV, we take as the allowed band of Ch a definition analogous to that used in 
Table I. This gives us the uncertainties in the corresponding predictions of thermodynamic quantities. As one could 
guess, the pleasantly stable results for thermodynamics shown in Figure 3 are model artifacts. In Figure 4 we show 
that the uncertainty bands in the two models defined by cqs. (3, 4) exclude each other in the vicinity of T = 170 McV. 
In view of these uncertainties, we do not consider the more detailed models which have been used [22]. In future, 
when observations of the hadron spectrum are extended to significantly larger masses, a more detailed consideration 
of string models may become a fruitful topic of research. 

Lattice computations of the equation of state give E/T^ = 1.4 ± 0.2 at a temperature of 140 MeV [15]. At such 
low temperatures the known hadron spectrum dominates the results (see Figure 4). The mismatch between lattice 
and hadron gas models is due to the fact that current day lattice thermodynamics computations are performed on 
lattices which are too coarse at such low temperatures. The best lattices have cutoff of around 1100 MeV when 
doing thermodynamics at T = 140 MeV. Since the cutoff is comparable to the low-lying baryon masses, the hadron 
resonance spectrum is strongly disturbed and the low-temperature results are not yet physical[37]. Since the lattice 
spacing varies inversely with the temperature, lattice results at higher temperatures are expected to be reasonable. 
The efficacy of the lattice closer to Tc is borne out by the fact that renormalization group invariant estimates of Tc 
are possible with the cutoffs in use today [17]. Interestingly, present day lattice computations at a temperature of 
210 MeV give E/T'^ = 12.4 ± 0.3 [15], which is also below the prediction of HRG2. However, as shown in Figure 4, 
at this temperature the problem very likely lies with the hadron gas models. We demonstrate this next. 

A phase transition involves a singularity in the free energy. As a result one usually does not have the same elementary 
excitations in terms of which the thermodynamics is constructed on both sides of a phase transition. However, there 
is no singularity of the free energy at a cross over. As a result, there is no theoretical bar to a model which is valid 
on both sides of the cross over. One way to investigate the applicability of such models on both sides of Tc is to 
look at quantum number susceptibilities (QNS). Introduce chemical potentials for the three quantities conserved in 
strong interactions, namely the baryon number, the third component of the isospin (J3) and the strangeness. Then 
the baryon number susceptibility is defined to be 



Xb{T) = 



d'^P{T, ij,B,i^i3,i^s) 



(5) 



This defines the susceptibility at zero chemical potential, since there is lattice data on this quantity, although, of 
course, it may also be taken at finite chemical potential. Following [18] the non-linear susceptibilities (NLS) are 
defined as 



d"P{T,HB,ni„fis) 



dill 



(6) 



The QNS is just the NLS for n — 2. In exact analogy one can also define strangeness susceptibilities and so-called 
off-diagonal susceptibilities where some of the derivatives are with respect to one chemical potential and others with 
respect to another [19]. The QNS predicted by HRGl and HRG2 are shown in Figure 5. The enhanced density of 
baryon states in HRG2 leads to a substantial increase in xb{T). As expected, current lattice data [24, 25] does not 
match these curves. 

The reason is very accurately probed by the quantity 



Cbs = -3 



XBsjT) 
XsiT) ■' 



(7) 



where xbs{T) is one of the off-diagonal QNS [26, 27]. The normalization is such that Cbs = 1 in an ideal gas of 
quarks. The HRGl and HRG2 predictions of Cbs are shown in Figure 6. Since the lowest mass baryons are non- 
strange, Cbs must start at zero. It is seen to climb monotonically with temperature. Cbs can exceed unity if the 
contribution of the singly and doubly strange baryons to xb exceeds half the contribution of strange mesons to xs- 
Since the strange meson spectrum starts at a much lower mass than the strange baryon sector, this cannot happen 
at low temperature. However, as seen in Figure 2, the observed density of strange baryons grows much more rapidly 
with temperature than that of strange mesons. As a result, at sufffciently high temperature Cbs exceeds unity. In 
HRG2 Cbs remains above unity until Th- Hence there is no continuity between this model and the physics of the 
high temperature phase of QCD. 
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FIG. 5: The baryon number susceptibility computed in hadron resonance gases. The first panel gives the results in HRGl and 
the second in HRG2. The colour coding is the same as in previous figures. Note the importance of including the Hagedorn 
spectrum of baryons. 
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FIG. 6: Cbs in the hadron resonance gas models HRGl (left) and HRG2 (right). One expects Cbs to be a monotonic function 
of temperature in a resonance gas model. 



The ratios of successive NLS in the HRG2 is shovifn in Figure 7. These ratios are not very sensitive to the change 
from HRGl to HRG2. Note that they are extremely constant as a function of temperature. In contrast, lattice 
computations show much structure in the vicinity of Tc as a consequence of the nearness of the QCD critical point 
[24]. The hadron resonance gas, being a mixture of ideal gases, sees no critical point, but only the Hagedorn limiting 
temperature. As a result, it misses all this structure. 

We have shown here that there are major points of mismatch between lattice computations of QNS and the resonance 
gas models for T > T^. Since there is no phase transition in QCD at hb =0 but only a crossover, the failure of a 
model soon above Tc also implies its failure a little before Tc- This means that the resonance gas models are restricted 
in their range of applicability to well below T^. In this real QCD seems to be very different from either pure gauge 
theory or large- iVc QCD, where the resonance gas model could remain perfectly accurate right up to the first order 
phase transition in these models. We turn next to the question of whether heavy-ion observables can give information 




FIG. 7: The non-linear baryon number susceptibilities computed in hadron resonance gas model HRG2. We have plotted the 
ratios x^"+'^^T^ /x^"^ for n = 2, 4, 6 and 8. The colour coding is the same as in previous figures. The results are very stable 
against the mass cutoff. 
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FIG. 8: Bjorken flow in the hadron resonance gas phase starting from an initial energy density of 1 GeV/fm^, with the bands 
of uncertainty corresponding to A = 1600 MeV, 1800 MeV and 2000 MeV for Cg. Also shown (in green) is the How for constant 
= 1/5. There is very little difference between these flows. 

on the Hagedorn spectrum of resonances. 

One of the uses to which equations of state can be put is hydrodynamics. We perform the following simple 
computation with the HRG2 equation of state — use it to evolve a longitudinally expanding fireball which has cooled 
to a point where its energy density, E, is 1 GeV/fm^. For longitudinal flow, the quantity E/Eq is a function only of 
the ratio t/tq where r is the proper time, and tq the initial time, when E{tq) = Eq. We integrated the longitudinal 
flow equation with upper and lower band of for A = 1600 MeV, 1800 MeV and 2000 MeV. We found that the effect 
of these changes is minimal, as we show in Figure 8. For comparison we also show the flow obtained with constant 
Cj = 1/5. This is almost indistinguishable from the other flows. Similarly, the result of a computation with the 
density of states in eq. (4) is indistinguishable from these. So, hydrodynamics is almost blind to the level of detail in 
the equation of state that we have studied. 

The main phenomenological application of hadron gas models, however, is in the analysis of particle yields in 
heavy-ion collisions. Most of the particle species observed in the detector are those which are stable under strong 
interactions, because the others decay long before reaching the detectors. The computation of the yield of particles 
in a hadron gas model is matched to data to extract the reaction volume and freezeout temperature and chemical 
potentials [5-7]. 

In HRGl one creates a table of decays from the particle data book. Each hadron Hi has decay modes labeled by 
a. The reaction 

j 

proceeds with a branching fraction B" , where nfj is the number of hadrons Hj produced in this reaction. For later 
convenience we add the trivial rule that stable particles decay to themselves with branching ratio unity. The expected 
number of Hj produced per decay of Hi are 

iV,,=^i?r<r (9) 

a 

For a stable particle Hi, we have Na — 1 and Nji = for all other j. In general, some of the decay products, Hj, 
may be unstable under strong interactions; in that case one has to follow the decay chain until only stable particles 
remain. The expected number of Hj resulting from Hi after all this, Afji, is easily found by sufficient number of 
matrix multiplications — 

J2 NJk,Nk,M■■■Nk,„,^. *-e- AA^iV™. (10) 

ki,k2,- --ykjn 

The minimum power m to which the matrix TV has to be raised is equal to the maximum number of steps in any decay 
chain. From the formula above it is clear that the rows of the matrix Af which correspond to unstable hadrons are 
zero. This means that Af = AfN, and hence if m is chosen to be larger than actually required, the cost is in CPU time 
and not in correctness. Proofs of all these assertions can be written down most simply by noting the structure of the 
matrix N, but insight is also gained by examining a one-to-one correspondence between this problem and problems 
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FIG. 9: The K/-k ratio obtained in our models using as input the energy dependence of the freezeout values of T and ^ obtained 
in [34] by fitting to the data exhibited. The first figure is for decay model 1 and the second for decay model 2. Only the results 
for A = 2000 MeV are shown for clarity. 



on directed graphs [28]. Finally, the expected yield of a stable particle, iJj, resulting from a fireball which freezes out 
at temperature T and chemical potentials ii is 

nii <A 

iV,(T/,T,A*)= ^ AG-.iVf(F,r,M), (11) 

i 

where Nf^ is the number of Hi in thermal equilibrium. From this, it is clear that for hadron yield computations, it 
not necessary to keep track of the detailed decays of every Hi. It is sufficient to keep A/i^ , i.e., the expected number 
of hadrons iJj, stable under strong interactions, resulting from the decay of of each Hi. In any case, the data may be 
read off the particle data book. 

In IIRG2 the yield computation is not so straightforward, since one has to start with a model Nji for all Hi. In this 
first study we restrict ourselves to models where Hj = tt, K or a generic baryon (i.e., we lump the ground-state octet 
and decuplet baryons into one). As a result, we cannot answer detailed questions about the yield, but only questions 
of total particle multiplicity and K /t: or mcson/baryon ratios. In general, we expect Mji (for fixed Hj) to depend on 
the mass of so, a generic model for decays is 

Nji — ai + a\mi + a^m^ + • • • . (12) 

While Mji may well depend on other quantum number of Hi, we are neither able to confirm or rule this out. Therefore, 
we decided to work with the simple model above in this paper. The decays of the heaviest hadrons are not yet studied 
well enough to provide further constraints. Further, incomplete data may contain unknown biases. Hence we neglected 
the heaviest hadrons and constrained the model with data for hadrons with masses up to 2000 MeV. 

Due to the high threshold for the production of a baryon-antibaryon pair, the decays of such particles do not involve 
such a pair in the final state. Such pairs are rare also in the known decays of hadrons with mass greater than 2000 
MeV, being seen mainly in the final state of strange mesons. The branching ratios to such pairs are still mostly ill- 
measured. As a result, we are unable to make any testable model about the expected occurrence of baryon-antibaryon 
pairs in decays of the Hagedorn spectrum of resonances. We therefore take the simplest possible model: 

^^B^ - B„ (13) 

that the number of baryons in the final state of the decay of any particle is equal to the baryon number of the 
resonance, Bi. When better data on resonances between 2000 and 3000 GeV becomes available, the model can be 
refined by adding higher terms from eq. (12). 

The production of KK pairs is also strongly suppressed in the decays of hadrons; this is a statement of the OZI 
rule. Cases which contradict this rule are known, but our attempts to make statistical models of such decays falls on 
the rock of insufficient data. Thus, we are forced to the rule that 

^fK^ = S,, (14) 

the number of kaons produced in the decay of Hi is equal to its strangeness, Si. We also examined a variant of this 
rule, which is that the strangeness in the baryon sector percolates down to strange baryons, and does not go into 
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production of kaons. This corresponds to the variant model Nki = Si{l — We refer to this as the decay model 

2, to distinguish it from the model in eq. (14), which we call decay model 1. 
Finally, we consider the expected number of pions produced in the decays of resonance. We fit the linear model 



where A^i is that part of the mass of Hi which is available for decay to pions. For unflavoured mesons, the available 
mass is the mass of the particle, for strange mesons, the available mass is the difference between the mass of a kaon 
and the mass of Hi, and for a baryon it is the difference between the mass of Hi and a typical baryon mass, which 
we take to be 1000 MeV. This simple model is fitted to four sets of data separately for strange and non-strange 
mesons and baryons. In all these cases we found h consistent with zero and a ~ 1.5 GeV~^, and consistent with each 
other. Adding in a term quadratic in masses to the model does not improve the fit significantly, so we keep to the 
linear model above. It may be useful to note the following implication: out of every GeV of rest mass of the higher 
resonances, about 500 MeV is available as the kinetic energy of each pion. 

This completes the specification of a simple model for the Mji in HRG2. Many more refinements and elaborations 
of the model are possible, however, this model is sufficient for the computations that we exhibit next. We extend the 
yield formula of eq. (11) to 



where wc use the empirically determined A/ij for < A and the model for mi > A. There are a very small number of 
baryons whose masses are sufficiently well-known to be included in HRGl, but whose decays are not very well known. 
For these we used the decay models. 

The K^/t:^ ratios have been determined at AGS [29], SPS [30] and RHIC [31] energies. In order to compare with 
this data we need to take into account the fact that the initial state has no strangeness. Due to strangeness conservation 
in strong interactions this condition of zero overall strangeness has to be enforced as a canonical constraint, i.e., the 
K/n ratios have to be determined in the canonical ensemble [32, 33]. We implement this in our computations of 
hadron yields. 

Our results for the K/ti ratio as a function of the beam energy \fS arc shown in Figure 9. In this computation we 
have used the freezeout parameters deduced in [34]. As expected, the ratios of yields are a little higher in decay model 
1. Note also that there is a little difference between the predictions of HRGl and HRG2. For the same freezeout 
conditions, IIRG2 gives a slightly smaller value of the yield ratio. In principle, the freezeout conditions that one 
deduces from data could be dependent on specifics of the hadron gas model one uses. Our results indicate that this 
dependence is at best mild. To check this we have implemented the freezeout criterion of [35] . Using this we find that 
the freezeout parameters change only marginally from HRGl to HRG2, consistent both with the conclusions of [36] 
and the results in Figure 9. In future we plan a more detailed study of hadron yields, and the effects of the Hagedorn 
spectrum on the chemical composition at freezeout in baryon rich matter. 

We conclude with a summary of our investigation into hadron resonance gas models of the kind which have been 
used to explain hadron yields in heavy-ion collisions. The thermodynamics of such resonance gases (HRGl) is strongly 
dependent on the cutoff in the spectrum at temperatures of over 100 MeV. In this paper we have investigated an 
extension of the hadron resonance gas models in which as yet unobserved resonances are included using the Hagedorn 
model of the hadron spectrum (HRG2) with two different densities of state (eqs. 3 and 4). We found that the Hagedorn 
temperature, Th, is strongly dependent on details of the model for the hadron spectrum. This implies that the present 
knowledge of the spectrum of QCD is as yet unable to constrain string models of hadrons. 

Each of the HRG2 models stabilizes the uncertainty in thermodynamics below due to the cutoff in HRGl. 
However, the predictions of the two variants of the model are different. While the model is defined up to T/f , we 
found that HRG2 gives unrealistic results for thermodynamics, especially for Cbs, from Tc to Th- Since there is 
no phase transition in QCD, a failure of a model near and above Tc also implies a failure near and below Tc- This 
leads us to believe that the accuracy of resonance gas models cannot be pushed much closer to Tc than the freezeout 
temperature. Of course, this leaves open the interesting possibility that resonance gas models describe the physics of 
transport that leads to freezeout, for example, the tliic;kness of the freezeout layer. 

We also investigated observable quantities in heavy-ion collisions, such as hadron yields. In order to do that, we 
had to develop a novel model for the decay of Hagedorn resonances. We found that the significant uncertainties that 
we saw in thermodynamics lead to rather small changes in observables such as the KJ-k yield ratio. This has positive 



f^-Ki = b + aMi 



(15) 




(16) 
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implications on efforts to explain hadron yields from heavy-ion observation. 
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